Finite Size Scaling of Domain Chaos 
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Numerical studies of the domain chaos state in a model of rotating Rayleigh-Benard convection 
suggest that finite size effects may account for the discrepancy between experimentally measured 
values of the correlation length and the predicted divergence near onset. 
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PACS numbers: 47.25.-c, 05.45.+b, 47.20.Tg 

Spatiotemporal chaos is the name given to states in 
driven nonequilibrium systems that are disordered in 
space and show chaotic time dynamics. The usual di- 
agnostic tools developed for chaotic dynamics in nonlin- 
ear systems with few dynamical degrees of freedom focus 
on the geometrical structures in the phase space of the 
dynamics. These largely lose their appeal for the very 
high dimensional dynamical systems of spatiotemporal 
chaos. The questions of how to characterize, understand, 
and predict the properties of spatiotemporally chaotic 
systems remain poorly understood, despite much exper- 
imental and theoretical attention over the past decade. 
Given this lack of understanding it is important to study 
systems which are favorable for both experimental and 
theoretical study, and in particular ones where theoret- 
ical predictions can be made and tested experimentally. 
In this paper we present results that further the compar- 
ison between theory and experiment for the state known 
as domain chaos |l| in rotating Rayleigh-Benard convec- 
tion. 

Rayleigh-Benard convection has long served as a 
canonical example of pattern formation in systems far 
from equilibrium. A fluid confined between two hori- 
zontal plates and driven thermally by maintaining the 
bottom plate at a higher temperature than the top plate 
undergoes an instability to a state in which there is fluid 
motion driven by the buoyancy forces induced by the 
thermal expansion. Far away from lateral boundaries 
the structure of the fluid motion locally forms the famil- 
iar convection rolls with a diameter close to the depth 
of the cell. This spontaneous formation of spatial struc- 
ture in a uniformly driven system is known as pattern 
formation. 

The instability to the roll state occurs at a particular 
value the Rayleigh number (a dimensionless measure of 
the temperature difference across the fluid) R — R c . For 
values of R slightly above R c the pattern is usually found 
to be a time independent state. However if the convection 
apparatus is now rotated about a vertical axis, Coriolis 
effects perturb the fluid velocity, and above a critical ro- 
tation rate Q c an ideal pattern of straight convection rolls 
is predicted to become unstable via the "Kuppers-Lorz" 
instability jl . The nature of the instability is that the 
rolls become unstable towards the growth of a second set 
of rolls with their axis rotated by an angle Okl from the 



axis of the original set of rolls. The value 9 k l depends 
on fluid properties, but for typical fluids is close to 60° at 
the onset of the instability. Once the second set of rolls 
grows to saturation replacing the first set, they in turn 
will become unstable towards rolls rotated through a fur- 
ther 9k Li etc. It is predicted that there will be no time 
independent saturated state even arbitrarily close to on- 
set. Experimentally the state of domain chaos is found, 
consisting of domains of differently oriented rolls with a 
persistent dynamics of domains expanding at the expense 
of one (or more) of the neighboring domains through the 
motion of the domain walls between them. 

The domain chaos state is of particular interest since it 
survives down to the onset of the convective state, where 
theoretical treatment based on an expansion in the weak 
nonlinearity should be possible. Indeed, based on the 
approximation that the switching angle is exactly 60°, 
Cross and Tu || developed a simple model of the domain 
chaos state, extending earlier work of Busse and Heikes 
||. The model is for the coupled dynamics of domains 
of rolls at only three orientations, characterized by three 
amplitudes Ai(x,y,t), i = 1,3 giving the strengths of 
the three roll components at each point in space and 
time. The coupled amplitude equations they used (after 
appropriate rescaling of space and time coordinates to 
eliminate unimportant constants) take the form 



d t A x = sA x 
d t A 2 = eA 2 
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Here the variable Xi is the spatial coordinate in the direc- 
tion perpendicular to the zth set of rolls; the derivatives 
transverse to these directions are of higher order. Al- 
though in general complex amplitudes should be used, 
Cross and Tu made the simplification of assuming the Ai 
to be real. This corresponds to neglecting variations of 
the wave numbers of the rolls. The parameter e measures 
the distance from onset (e = (R — R c )/R c ), and is the 
small parameter of the expansion. The constants g + and 
<?_ determine the interaction between one set of rolls and 
the set rotated through +60° and —60° respectively. In 
the absence of rotation, clockwise and anticlockwise rota- 
tions are equivalent so that g + = g_ and in this limit it is 
easily shown that Eqs.(P are relaxational or "potential" 
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H so that persistent dynamics is impossible. The rota- 
tion breaks this "chiral symmetry" , so that g + ^ g_ , and 
the equations are then no longer relaxational. Cross and 
Tu found numerically, for sufficiently different g + , <?_ , a 
state of persistently dynamic domains. 

The demonstration of the existence of the chaotic do- 
main state within Eq. (|l|) gives immediate quantitative 
predictions for the scaling of the size and lifetime of the 
domains with s. Defining scaled variables Xi = e 1 / 2 :^, 
Ti = et, and Ai = e~ 1 / 2 At, yields equations with no ap- 
pearance of the parameter e. In these scaled variables 
the domain size and lifetime are therefore e independent, 
leading to the prediction of a domain size or correlation 
length £ of the state scaling as e -1 / 2 and a domain life- 
time or correlation time r of the dynamics scaling as e _1 
in the physical variables. This remains one of the few 
quantitative predictions for properties of spatiotcmpo- 
ral chaos in an experimental dissipative system far from 
equilibrium that has actually been tested experimentally. 
However recent experiments || did not find the predicted 
result. It is this important discrepancy that we address 
in the present paper. 

The clear disagreement between the experiment and 
the theory is that in the experiment £ -2 and the domain 
switching rate u> a = r _1 appear to remain finite as e ap- 
proaches zero, instead of going to zero linearly in e. If a 
power law dependence on s is forced on the data, powers 
much smaller than 1 result. Based on numerical simula- 
tions of equations showing domain chaos we suggest that 
this discrepancy between experiment and theory might 
be due to the (necessarily) finite size of the experimental 
system. 

Our conclusion is based on the following. We find re- 
sults for the correlation length in numerical simulations 
that have many of the same qualitative features as in 
the experiment. In particular the measured correlation 
length £m (using an algorithm defined below that is the 
same as the one used in the experimental work) appears 
to remain finite approaching the threshold. The flexi- 
bility of the numerical approach allows us to investigate 
this behavior as a function of the aspect ratio T. We find 
that the data for different T and different e collapse onto 
a single form suggested by finite size scaling 

6/ - e/(e/r) (2) 

with £ oc e -1 / 2 the "ideal" correlation length following 
the theoretical prediction. Our numerical data is con- 
sistent with f(x) — > const for small x so that in large 
enough systems (or for large enough e) the predicted 
dependence proportional to £ ~ e -1 / 2 would be found, 
whereas f(x) oc x^ 1 for large x, so that in small sys- 
tems (or for small e) the measured correlation length is 
proportional to the system size. 

The equation we simulate is the partial differential 
equation for a real field tp(x,y,t) in a two dimensional 



domain 

9 t V> = e^) + (V 2 + I) 2 V> - SiV* 3 

+ g 2 z- V x [( W) 2 W] + 33 V- [(V^) 2 W] • (3) 

with boundary conditions 

ip = n • W = (4) 

on the boundary with normal n. This equation has been 
introduced previously to model domain chaos ]7j] with 
periodic boundary conditions. With g 2 = 33 = it 
is the well known Swift-Hohcnberg equation that has 
been much studied in the context of pattern formation 
||. The spatially uniform state -0 = becomes unsta- 
ble to a stripe state with wave number q c = 1 at e = 0, 
and for e > steady, stable, nonlinearly saturated stripe 
states may be found. In modelling a convection system 
we can imagine tp(x,y) as representing the temperature 
field across a midplane of the system, and the stripes 
are a section of the convection rolls. The second nonlin- 
ear term in Eq. (^|), with coefficient gi is the all impor- 
tant term yielding the breaking of the chiral symmetry 
induced by rotation in the physical system — and so in- 
creasing 172 corresponds to increasing rotation rate in the 
convection system. For sufficiently large gi the stripe 
state becomes unstable to a rotated set of stripes as in 
the Kuppers-Lorz instability. The angle at which the new 
set of stripes occurs can be tuned with the parameter g% 

Although Eq. mf) is easier to evolve numerically than 
the coupled equations for fluid velocity and temperature 
fields in a three dimensional domain that give a complete 
description of the convection system, the task remains 
challenging. This is because we need to integrate the 
equation in a large domain and over long times. In fact, 
since we are interested in approaching the limit e — > to 
uncover the scaling behavior, and, in this limit, the dy- 
namics becomes very slow, exceedingly long integration 
times are necessary. In addition, to model the experi- 
ment the domain must be circular, and to eliminate any 
bias towards a particular stripe orientation the use of cir- 
cular polar coordinates to describe the geometry seems 
preferable. 

To meet these numerical challenges we have developed 
a fully implicit method using a finite difference represen- 
tation of the differential operators on a polar coordinate 
mesh. At each time step the nonlinear Crank-Nicolson 
equations for the new value of the field are solved by 
Newton's method. The GMRES (generalized minimal 
residual) iterative method preconditioned with the 
Bjorstad || fast direct biharmonic solver, is used to solve 
the linear problem within each Newton iteration. This 
method allowed large time steps Ai, for example up to 
100 at e = 0.01. A variable time stepping algorithm 
was implemented to exploit this opportunity, based on a 
comparison of results with time steps At and At/2. This 
time stepping should be compared with what might be 
obtained with a conventional semi-implict method where 
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the gradient terms in the nonlinear terms impose a sta- 
bility limit on the possible time step determined by the 
spatial mesh resolution. This limitation is particularly 
severe using polar coordinates since the mesh spacing in 
the azimuthal direction becomes very small near the po- 
lar origin. In practice we found a time step limited to 
At < 0.1 at e = 0.01 using a conventional semi-implicit 
code, a factor of 1000 smaller than in the fully implicit 
code. 

The polar mesh induces an artificial singularity in the 
description at the origin of the circle where the physical 
behavior will be smooth. This difficulty has plagued the 
development of many codes based on polar coordinates. 
To test the accuracy of our code in the vicinity of the ori- 
gin we simulated Eq. (Q) with gi = g 2 = starting with 
an initial condition of straight stripes. With these pa- 
rameters the stripe state is unstable to a square pattern. 
For an initial condition that is exactly a stripe state, the 
development of squares can physically only be initiated 
at the boundary, and the square pattern should then be 
observed to propagate in from the boundaries. With in- 
adequate spatial resolution we found instead that squares 
also began to grow via nucleation around the origin, pre- 
sumably an artifact of the reduced accuracy of the nu- 
merical integration here. This unphysical result could be 
effectively eliminated by increasing the number of radial 
mesh points (as mentioned above, the azimuthal resolu- 
tion for a fixed number of azimuthal points becomes very 
fine near the origin) . To take the best advantage of the 
computer resources we used a variable mesh in the radial 
direction, with the resolution increasing smoothly to a 
factor of two or three improvement over the inner third 
of the radial direction. It should be noted that these test 
runs use the exponential amplification due to a physical 
instability (stripes to squares) acting over a long time 
(e.g. a time of 100) to enhance the effects of numerical 
inaccuracies to visible amplitudes. Thus the elimination 
of visible spurious effects near the origin in these test 
runs gives us confidence in the accuracy of the numerical 
procedure. 

To study domain chaos we simulated Eq. (||) with fixed 
values of the nonlinear coefficients g\ = 1, g% = —2.60, 
and 53 = 1.5 corresponding to a Kuppers-Lorz instability 
at an angle of 60° |7]] . We studied the behavior for values 
of £ between 0.01 and 0.3 and in circular geometries with 
radii 307r, 407r, 50-7T, and 807r, corresponding to aspect 
ratios of 30 — 80 in the experimental system. Snapshots 
of the field tf)(x,y) from two runs are shown in Fig. (Q). 
Notice the dependence of the domain size on the control 
parameter e, and the comparable values of domain and 
system sizes at the smaller value of e. 

The quantitative diagnostics of the state are based on 
the Fourier transform ip()i) using a Hanning window over 
an inscribed square, following the same approach as in 
the experiments ]1(|. The intensity S(k) = |-0(k)| 2 is 
found to be concentrated on a ring in Fourier space near 
k = 1 corresponding to the stripe periodicity of 2-7T. The 
intensity varies with the angle around this ring, corre- 




FIG. 1: Snapshots of the dynamic domain configuration at 
two values of the control parameter (a) e = 0.01 and (b) 
e = 0.3. The aspect ratio is 40n. 



sponding to the varying representation of domains of dif- 
ferent stripe orientations in the pattern. This angular 
distribution varies with time as the domains evolve, and 
can be used to define the switching rate uj a . The correla- 
tion length of the pattern is defined as the average width 
of the ring in Fourier space £m = [(k 2 ) — (k) 2 ]^ 1 / 2 with 



Jk n S(k)d 2 k 
JS(k)d 2 k 



(5) 



where () t denotes a time average. 

Results for the correlation length are shown in Fig. <^). 
As in the experiments the inverse of the measured corre- 
lation length £Tj appears to remain nonzero at onset, and 
the data are reasonably well fit by a form ~ y/e — So 
rather than the expected yfe dependence. The depen- 
dence of Eq wc find on the aspect ratio suggests that 
this might be understood as a finite size effect. This 
is confirmed by the scaling plot motivated by Eq. (|^) 
which suggests that a plot of e -1 / 2 /^/ against e~ 1 ^ 2 /T 
should collapse the data. The successful collapse is shown 
in Fig. (||). Note that for small values of e^ 1 / 2 /T (i.e. 
large T or e) the value of ^Af£ -1 ^ 2 approaches a constant 
corresponding to the theoretical expectation for large 
enough systems. On the other hand for large e~ 1 ^ 2 /T the 
curve approached a linear behavior consistent with the 
correlation length scaling simply with the aspect ratio, 
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FIG. 2: Inverse measured correlation length f^ 1 as a function 
of control parameter e for various aspect ratios F. Fhe lines 
are fits of the form = a(e — eo) 1//2 - 




k — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — L 
0.05 0.10 



FIG. 3: Plot of the inverse correlation length £ A / scaled with 
£ 1//2 against the inverse aspect ratio F _1 scaled by e -1 ^ 2 
showing the collapse of data for different e and F onto a 
single curve. Fhe solid line is a fit to the empirical form 
y = Va 2 + b 2 x 2 yielding a = 0.091 and b = 2.8. The dashed 
lines show the asymptotic limits of this fit. 



which is much harder to do experimentally than numer- 
ically. The experiments did investigate the dependence 
of the correlation length on the rotation rate, and found 
that the measured correlation lengths at different rota- 
tion rates could be related by an e-independent scale fac- 
tor: £(e,Q) ~ £o(n)f(e). Since the deviation from the 
expected e" 1 / 2 behavior found in /(e) is common to all 
the runs at different fi, this scaling does not shed light 
on the basic discrepancy with the predicted e" 1 / 2 scaling 
however. 

We have also studied the behavior of the switching 
frequency to a on e. Here, unlike the experiment, we find 
good agreement with the prediction oj a oc e for small 
e with no dependence on the aspect ratio. The trend 
c<j a — ^ for e — > even in finite size systems is not surpris- 
ing for Eqs.(||[I]) since the effects of the rotation that lead 
to the persistent dynamics only appear in the nonlinear 
terms which become small for e — > 0. A careful analy- 
sis of the fluid equations shows that there are also linear 
terms depending on the rotation that are important near 
the boundaries. This means that in a finite system the 
onset of convection rolls in a rotating system is actually 
a Hopf bifurcation with the onset frequency going to zero 
for large aspect ratios, and so it might not be surprising 
that co a remains nonzero at onset. These effects can in 
fact be captured using modified boundary conditions for 
the Swift-Hohenberg equations . However these com- 
plicated boundary conditions make the numerical scheme 
considerably harder, and we have not modified our code 
to investigate them. 

In conclusion our numerical simulations of model equa- 
tions for rotating Rayleigh-Benard convection show re- 
sults for the correlation length near threshold that show 
qualitative similarities to the experimental results. A fi- 
nite size scaling ansatz shows that our measured lengths 
are in fact consistent with the expected e -1 / 2 divergence 
near threshold, but that the finite system size obscures 
this dependence. It will be interesting to see if this re- 
sult can be confirmed experimentally, or in simulations 
of the full fluid and heat equations for convection that 
we arc currently pursuing. It is noteworthy that similar 
predictions for a diverging correlation time near thresh- 
old in electroconvection spatiotemporal chaos have re- 
cently been verified experimentally |l2| . In this system 
larger aspect ratios are accessible so that finite size effects 
should not be important. 



£m — r/2.8. Note that the finite size corrections become 
important (e.g. identified as the intersection point of the 
straight lines in Fig. (|h for T < 3£, with an additional 
factor of 3 over the naive expectation. 

In the experimental work there has been no system- 
atic attempt to study the dependence on aspect ratio, 
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